Load packages, functions, paths

library(dplyr)
library(tidyr)
library(ggplot2)
library(calibrate)
library(knitr)
"%&%" = function(a,b) paste(a,b,sep="")
my.dir = "/Volumes/dolan-lab/hwheeler/ThePlatinumStudy/GWAS/"
res.dir <- my.dir %&% "genotypes/UMich_imputation_results/PrediXcan/PrediXcan_results/"
source(my.dir %&% "genotypes/qqunif_px.r")
source(my.dir %&% "manhattan.R")

ord3CIPN8 = pred_exp + agediagnosis

Bonf. sig. tissues:

  • DGN-WB
  • TW_Brain-Nucleusaccumbens-basalganglia_ElasticNet.0.5
  • TW_Breast-MammaryTissue_ElasticNet.0.5
  • TW_Cells-Transformedfibroblasts_ElasticNet.0.5
  • TW_SmallIntestine-TerminalIleum_ElasticNet.0.5
  • TW_Spleen_ElasticNet.0.5
  • TW_Stomach_ElasticNet.0.5
  • TW_WholeBlood_ElasticNet.0.5
pheno <- 'ord3CIPN8'
tislist <- scan(my.dir %&% "pred_exp_list","c")
for(tis in tislist){
  res = read.table(res.dir %&% "N88_n953_ord3CIPN8_agediagnosis_" %&% tis %&% ".ordreg.PrediXcan",header=T)
  newres = mutate(res,CHR=as.numeric(gsub("chr","",chr)),BP=start,SNP=gene) %>% arrange(CHR,BP)
  res <- newres[complete.cases(newres),] #rm NAs
  res <- res[res$P != 0,]#rm P=0 (unstable, see permutations below)
  if(min(res$P)<0.0001){
    manhattan(res,suggestiveline=FALSE,genomewideline = -log10(0.05/dim(res)[1]),main=pheno %&% " " %&% tis,annotatePval = 0.0001,cex.main=0.8)
    tophits <- dplyr::filter(res,P<=0.0001) %>% arrange(P) %>% dplyr::select(-CHR,-BP,-SNP,-lCI, -uCI)
    cat('\n')
    print(kable(tophits))
    cat('\n')
  }
  else{
    manhattan(res,suggestiveline=FALSE,genomewideline = -log10(0.05/dim(res)[1]),main=pheno %&% " " %&% tis,cex.main=0.8)
    cat('\n')
  }
  qq<-qqunif(res$P,plot=TRUE,title=pheno %&% " " %&% tis)
  cat('\n')
}

Value Std..Error t.value OR n gene P chr strand start end ensid
-2.772388 0.5954898 -4.655642 6.251260e-02 677 RPRD1B 3.20e-06 chr20 + 36661948 36720768 ENSG00000101413.7
14.273709 3.3613674 4.246399 1.581222e+06 677 GAB2 2.17e-05 chr11 - 77926343 78129394 ENSG00000033327.8

Value Std..Error t.value OR n ensid P chr strand start end gene
11.42955 2.890188 3.954604 92000.34 677 ENSG00000101347.7 7.67e-05 chr20 - 35526225 35580246 SAMHD1

Value Std..Error t.value OR n ensid P chr strand start end gene
-4.290414 1.030405 -4.163815 0.0136993 677 ENSG00000101413.7 3.13e-05 chr20 + 36661948 36720768 RPRD1B

Value Std..Error t.value OR n ensid P chr strand start end gene
-0.6445523 0.1655098 -3.894347 0.5248975 677 ENSG00000110455.9 9.85e-05 chr11 + 44087475 44105764 ACCS

Value Std..Error t.value OR n ensid P chr strand start end gene
-4.809874 1.075943 -4.470379 0.0081489 677 ENSG00000103653.12 7.8e-06 chr15 + 75074398 75095539 CSK

Value Std..Error t.value OR n ensid P chr strand start end gene
3.810047 0.8740653 4.358996 45.15256 677 ENSG00000174606.8 1.31e-05 chr1 - 213165524 213189168 ANGEL2

Value Std..Error t.value OR n ensid P chr strand start end gene
-3.520886 0.6956810 -5.061064 0.0295732 677 ENSG00000101413.7 4.00e-07 chr20 + 36661948 36720768 RPRD1B
-1.583702 0.3877443 -4.084397 0.2052140 677 ENSG00000186314.7 4.42e-05 chr5 - 144851362 145214899 PRELID2

Value Std..Error t.value OR n ensid P chr strand start end gene
-0.0081341 0.0009698 -8.387602 0.9918989 677 ENSG00000101544.4 0 chr18 + 77866915 77905406 ADNP2

Value Std..Error t.value OR n ensid P chr strand start end gene
1.405081 0.3352211 4.191505 4.075857 677 ENSG00000142192.16 2.77e-05 chr21 - 27252861 27543446 APP

Value Std..Error t.value OR n ensid P chr strand start end gene
-5.007081 1.150917 -4.350516 0.0066904 677 ENSG00000183914.10 1.36e-05 chr17 + 7620672 7737062 DNAH2

Value Std..Error t.value OR n ensid P chr strand start end gene
-3.453223 0.8585585 -4.022118 0.0316435 677 ENSG00000110680.8 5.77e-05 chr11 - 14988214 14993900 CALCA

Value Std..Error t.value OR n ensid P chr strand start end gene
-10.0862840 2.4318382 -4.147597 0.0000416 677 ENSG00000152127.4 3.36e-05 chr2 + 134877554 135212192 MGAT5
-1.6455843 0.4051854 -4.061312 0.1928998 677 ENSG00000101407.8 4.88e-05 chr20 - 36611409 36661870 TTI1
0.8512563 0.2105865 4.042311 2.3425880 677 ENSG00000132821.7 5.29e-05 chr20 + 36531499 36573752 VSTM2L
0.9169173 0.2320064 3.952120 2.5015670 677 ENSG00000176715.11 7.75e-05 chr16 + 89154783 89222254 ACSF3

Value Std..Error t.value OR n ensid P chr strand start end gene
14.41464 3.592181 4.012782 1820531 677 ENSG00000186480.8 6e-05 chr7 + 155089486 155101945 INSIG1

Value Std..Error t.value OR n ensid P chr strand start end gene
1.172545 0.2661875 4.404959 3.230202 677 ENSG00000126262.3 1.06e-05 chr19 + 35934809 35942667 FFAR2
4.273914 1.0501990 4.069623 71.802121 677 ENSG00000124882.3 4.71e-05 chr4 + 75230860 75254468 EREG

Value Std..Error t.value OR n ensid P chr strand start end gene
-2.936122 0.7536232 -3.896008 0.0530711 677 ENSG00000164877.14 9.78e-05 chr7 - 1473996 1499138 MICALL2

Value Std..Error t.value OR n ensid P chr strand start end gene
1.636136 0.3635916 4.499929 5.135291 677 ENSG00000057019.11 6.8e-06 chr3 - 98514785 98620533 DCBLD2

Value Std..Error t.value OR n ensid P chr strand start end gene
-4.857885 1.065442 -4.559501 0.0077669 677 ENSG00000101413.7 5.1e-06 chr20 + 36661948 36720768 RPRD1B

Value Std..Error t.value OR n ensid P chr strand start end gene
2.183234 0.4560312 4.787465 8.874957 677 ENSG00000132821.7 1.7e-06 chr20 + 36531499 36573752 VSTM2L

Value Std..Error t.value OR n ensid P chr strand start end gene
-3.55916 0.799912 -4.449439 0.0284627 677 ENSG00000177156.6 8.6e-06 chr11 + 747329 765024 TALDO1

PrediXcan 10,000 perms select tissues

tislist=c('DGN-WB','TW_Brain-Nucleusaccumbens-basalganglia_ElasticNet.0.5','TW_WholeBlood_ElasticNet.0.5')
tis="DGN-WB"
for(tis in tislist){
  overall<-read.table(res.dir %&% 'perms/N88_n953_ord3CIPN8_agediagnosis_' %&% tis %&% '.ordreg.PrediXcan.empP_permset1',header=T) %>% arrange(ensid)
  count.mat <- matrix(NA,nrow=dim(overall)[1],ncol=100)
  for(i in 1:100){
    res <- read.table(res.dir %&% 'perms/N88_n953_ord3CIPN8_agediagnosis_' %&% tis %&% '.ordreg.PrediXcan.empP_permset' %&% i,header=T) %>%   arrange(ensid)
    count.mat[,i]<-res$counta
  }
  empP <- rowSums(count.mat)/10000
  newoverall<-cbind(overall[,1:14],empP) %>% mutate(CHR=as.numeric(gsub("chr","",chr)),BP=start) %>% arrange(CHR,BP)
  print(ggplot(newoverall,aes(x=P,y=empP))+geom_point(alpha=1/5)+ggtitle(tis %&% ' PrediXcan 10k perms'))
  res <- newoverall %>% mutate(empP=ifelse(empP==0,1e-05,empP),SNP=gene) #replace empP=0 with 1e-05
  res <- res %>% dplyr::filter(P>0) #remove P=0 (unstable)
  res <- res[complete.cases(res),] #rm NAs
  kable(head(arrange(res,empP)))
  manhattan(res,suggestiveline=FALSE,genomewideline = -log10(0.05/dim(res)[1]),main="P ord3CIPN8 " %&% tis, annotatePval = 0.0001)
  qq<-qqunif(res$P,plot=TRUE,title="P ord3CIPN8 " %&% tis)
  empPres <- mutate(res,P=empP)
  manhattan(empPres,suggestiveline=FALSE,genomewideline = -log10(0.05/dim(empPres)[1]),main="empP ord3CIPN8 " %&% tis, annotatePval = 0.0001, ylim=c(0,5.2))
  qq<-qqunif(res$empP,plot=TRUE,title="empP ord3CIPN8 " %&% tis)
}